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Chapter 1 


INTRODUCTION 

The steady operation of a liquid-propellant rocket 
engine is often disturbed by the occurrence of large 
pressure oscillations in the combustion chamber. These 
oscillations, which can lead to damage to or failure of the 
motor, are caused by amplification of initially small 
acoustic disturbances due to the energy released by unsteady 
burning of the propellant. This is usually referred to as 
combustion instability. It is generally accepted that 
unsteady burning can be correlated with the gas pressure 
(pressure sensitivity) and the magnitude of velocity of the 
fuel drops relative to the gas (velocity sensitivity). 

This dissertation is concerned with mathematical modeling 
of velocity-sensitive combustion instability in liquid- 
propellant rocket motors. 

A survey of literature dealing with mathematical 
modeling of pressure-sensitive combustion instability can 
be found in the dissertation of Powell (1). One of the 
first papers to examine nonlinear effects is that of Maslen 
and Moore (2) who considered only the fluid mechanical 
effects in a circular cylinder. The paper by Priem and 
Guentert (3) is one of the first to include the effect of 
velocity sensitivity. They discussed combustion instability 
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in a thin annulus . 

Since direct numerical solution of the governing 
equations is very time consuming and difficult to extract 
information from, many Investigators such as Powell and 
Zinn (4), Lores and Zinn C5)j and Peddiesion, Ventrice, and 
Purdy (6) have employed methods of weighted residuals such 
as Galerkin and collocation methods . 

Previous applications of the method of weighted 
residuals to combustion-instability problems have dealt with 
pressure-sensitive combustion. In this work this method will 
be extended to handle situations in which velocity sensitivity 
is important. Both the collocation and Galerkin methods will 
be considered. Stability boundaries and pressure wave forms 
will be computed numerically for transverse motion in a - 
cylindrical combustion chamber. In addition a finite- 
difference method will be used to solve the one-dimensional 
problem of transverse motion in a thin annular chamber. 

These results will be compared with those obtained by a 
method of weighted residuals solution of the same problem 
in order to assess the accuracy of the approximate solution. 

In this study, attempts to solve the velocity- 
sensitive combustion instability problem by various 
numerical methods are considered. The governing equations 
that describe the flow conditions Inside liquid-propellant 
rocket motors will be derived in Chapter Two. An analytical 
solution is obtained in Chapter Three for nonlinear acoustic 
motion in a cylindrical chamber. The collocation method is 



applied to combustion instability in a cylindrical chamber 
in Chapter Four. In Chapter Five, the Galerkln method is 
applied to the same problem. In Chapter Six, attention is 
given to the analysis of combustion instability in a thin 
annulus. This allows the comparison of the Galerkln method 
and a finite-difference solution in a one-dimensioned 
context. Finally, a summary of this research is contained 
in Chapter Seven. 



Chapter 2 


COMBUSTOR EQUATIONS 


It is the objective of the discussion presented in 
this chapter to analyze the combustor conservation equations 
and reduce them to a tractable system. The simplification 
will be done in such a manner that will allow the resulting 
equations to retain both the mathematical and physical 
essence of the original problem. 

Under the assumptions that the usual balance 
principles of mechanics can be applied separately to each 
phase, the following equations are derived by applying the 
laws of conservation of mass, momentum, and energy to an 
arbitrary stationary control volume. 

Conservation of mass of the fluid requires that the 
time rate of change of mass in a volume v equals the rate 
at which mass enters v plus the rate at which mass is 
generated inside v and can be expressed as follows (for a 
fixed volume): 

* *->•->• * 
pdv ~ - f pu*nds + / wdv 

^ V S' V 


* 

where n is an outward unit vector normal to s , p is the 

^ * 

density of fluid, u is the velocity of fluid, w is the rate 
at which mass is generated per unit volume, and t is the time. 
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Using th.e divergence theorm and the arbitrariness of 
the stationary control volume, yields 

+ v»(pu) = w . (2.1) 

Similarly, the balance law of mass for fuel phase is 

* * 

+ V*(p^u^) = - w (2.2) 

* : 
where p is the density of fuel and u is the velocity of 
L/ L 

fuel . 

The balance law of linear momentum states that the 
rate of change of linear momentum in the volume v equals 
the rate of entry of linear momentum across the surface s 
plus the sum of all external forces , acting on the volume v 
and can be expressed as 

» „ 

3 a/ pudv = - / pu(u*n)ds + / Fdv + / wu-rdv - / Pnds 
t" V s V V ^ s 

* 

v/here is the force per unit volume applied to the gas by 
the fuel and P is the pressure of the gas. 

Using the same argument, the equation becomes 

3^*(pu) + ^*(puu) = ? + wu^ - VP. (2.3) 

Similarly, the balance of momentum for the fuel phase is 



6 


* »* 

- wu . (2.4) 

L 

Substituting C2.1) in (2.3) yields 

p(9^^u + u*^) = f + wCuj^ - u) - ^P. (2.5) 

Similarly for the fuel phase 

= - P. (2.6) 

The law of conservation of energy states that the 
rate of change of energy in a volume v equals the rate at 
which energy enters the volume v plus the rate at which 
energy is generated internally plus the rate at which work 
is done by external forces and gives 


t* 


» * 


* * * * 

= - 


3^jf/^p(e + ^u^)dv = 


* * * 

- / p (e + 35U^)(u'n)ds - / (Pn)*uds 
s s 


* ^ 
f 


V 


* .» 


u^dv + f w(ej^ + 3gu£)dv + / Qdv 


where e + hn^ is the energy of gas per unit volume, e + 5^u^ 

L/ 1j 

¥: 

is the energy, of fuel per unit volume, and Q is the heat 
transfer rate from the fuel to the gas. 


Following the same argument, the equation becomes 



7 


# ft „ «■ 

9 jfCpCe + h^^)) + v*(p(e + hn^)\x) = -V • (Pu) 


* » 




(2.6) 


Similarly i the balance of energy for the fuel phase becomes 


* * 


%U2 ) ) 
Li 


^ ^ 


* * 


* * 



* » 
w(e^ 


h^i) 



(2.7) 


Substituting (2.1)j (2.2) into (2.6), (2.7) respectively 
yields 


p (3 »(e + ) + 

U 


u 


■VCe 


+ 5su2)) 


V • (.ru; 


■u + w((e + ^u2) _ (e + %u2)) + Q 
ij Li Li 




( 2 . 8 ) 




* (©1 


hSp 


t-^cl 


+ ^U 2 ) ) = - i . U - Q . 
Li L 


(2.9) 


It is more convenient to work with the thermal energy 


equation for each phase. For gas phase this is done by 
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* 

dotting u into (2.3) to obtain tlie meclianical energy 
equation and subtracting this from (2.8) and using (2.1) to 
get 


* * * * * 
pDe/Dt*= PD /Dt*(logp) + ?*(Ut - u) + Q + w((e_ - e) 

Lj L 

+ (u^ - u )/2 + u*(u - u^) - P/p) . (2.10) 

In a similar way, it can be shown that 

* * * » , . 
p D e /Dt = - Q (2.11) 

L L L 

* * * 

where D /Dt* = 3 *+ u*?, D^ /Dt* = 3 + u^. are the 

t* L t* L 

comoving derivatives for each phase. 

For simplicity, gas phase viscosity and heat 
conduction were neglected in the previous analysis. Since 
these produce dissipation, equations derived in this way 
should give a conservative estimate of the stability of the 
system. 

It appears {see, for instance, Powell ( 1 )} that the 
primary effect of the burning fuel is that of interphase 
mass transfer. Thus, the balance laws for mass, momentum, 
and thermal energy for each phase will be further 
simplified by neglecting all Interphase transfer terms 
other than those appearing in the mass-balance equations. 
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This leads to the system of equations: 

* , * -V * 

Dp/Dt + pv*u = w 
* 

Dp /Dt + p_v*u_ = - w 

D Li D 

* *» 
pDu/Dt* = - 

* * * 
pLOhU^/Df = 0 


pC^DT/Dt* = P(D /Dt*)(logp) 

* * 

p C DT /Dt*= 0 (2.12) 

L vii L 

* » » •» 

where the constitutive equations e = C^T and e^ = ^vl'^L 
have been used (vjlth T denoting temperature and specific 
heat). The equation of state for an Ideal gas Is 

* 36 * 

P = pRT (2.13) 

where R Is the gas constant. 

The governing equations will now be nondlmenslonallzed 
with respect to a steady reference state, which will be 
denoted by the subscript "r” . All lengths will be referred 
to some characteristic length , such as the chamber radius. 
The characteristic velocity Is the speed of sound at the 



reference state, and the characteristic time Is the wave 


travel time The 

defined below; 

* 

V =v /L^ 

* 

u = C^u 

dimensionless quantities are 

! = Lj,t/C^ 

®L= 

3t 


P = PpP 

PL“ 

* 


P = PpP 

w = p^C^w/Lj, 

« 

* 

T = c^^t/CyR) 

2 

Tl= C^Tl/(yR) 

2 

2 

Cp= YPp/ Pp 

Pp“ PpP^"^! * 

The dimensionless conservation equations become 


Continuity 

Dp/Dt + pV*u = w 


Dpj^/Dt + pj^^*u^ = -w (2. 


Momentum 


Du/Dt = - ^P/y 
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Du^/Dt = 0 


Energy 


DT/Dt = (y - 1)P (D /Dt)(logp) 


Equation of state 


DT^/Dt = 0 


P = pT . 


By solving (2.16a) one obtains 


(2.15) 


(2.16) 


(2.17) 


T = pY-i . (2.18) 

Equation (2.17) then becomes 

P = P^ . (2.19) 

Substituting (2.19) to (2.15a) produces 


Du/Dt = -V(p^-1)/(y~ 1) . (2.20) 

Taking the curl on both sides of equation (2.20) gives 

^ X Du/Dt = 0 . (2.21) 

It can be shown using vector identity that (2.21) 
leads to the following equation which describes the 
generation of the vorticity 0 = ^ X u in the flow. 


Dfi/Dt = fi-(v u) - ^(V*u) . 


( 2 . 22 ) 
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This equation describes the variation of vorticlty 
for a given fluid particle. Suppose that at t = 0 j this 
particle has zero vorticlty and at the point where It Is 
located the velocity gradients are bounded. Then by (2.22) 
the rate of change of vorticlty of the particle vanishes. 

It follows from this that the vorticlty will be zero at the 
next Instant. By induction, it is seen that as long as the 
velocity gradients are bounded at each point occupied by the 
particle, its vorticlty will remain zero for all time. If 
all fluid particles in the system have zero vorticlty at 
t = 0, the vorticlty will vanish at all points in the flow 
field and for all t^O, therefore, as long as the velocity 
gradients are bounded through the flow field ( as they would 
not be, for example, at a shock wave ). Assuming this 
condition to be satisfied Implies Irrotatlonal flow (fi = 0) 
and allows the definition of a velocity potential (j) defined 
such that 


u = V<p . (2.23) 

Substituting (2.23) to (2.20) gives 


9 ^<{) + +pT~^/ ( y- 1) = F(t) . 

Since the velocity is the space derivative of (p , it 
is permissible to add to ^ any arbitrary function of time. 
This is equivalent to a statement that F(t) is arbitrary. 
For future convenience F(t) was selected to be 1/(y-1). 
This choice yields 
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pY-l= 1 _ Cy-1)C3^<J> + . (2.24) 

Substituting (2.24) into (2.l4a), (2.l8) and (2.19) gives 

3,,^ - V^^(l - (y~ 1)(3,<I' + ^V<()*V(j))) + V(j)*V(23 d) 
tt t t 

+ ^V^-V^) + (1 - (y- 1)(3 ((> + JgV(}.*V(|)) ) ^Y“2 )/(y-1)^= 0 

X/ 

T = 1 — (y— 1)(3 (|) + ^^(f)*^({)) 

P = (1 - (Y-l)(3^«j> + . (2.25) 

t 

This system of equations contains nonlinear gas 
dynamics terms of all orders. 

Due to their highly nonlinear and mathematically 
complicated nature, the system of equations obtained above 
cannot be solved exactly. In order to obtain simpler, but 
approximate, equations which can be more easily analyzed, all 
nonlinear terras of order higher than two are neglected. 

This has the effect of retaining the most Important 
nonlinear effect, while greatly simplifying the form of the 
equations. The system of equations then becomes 

p = 1 - (9^<i> + + ^(2 -y ) ( 9^(j> ) ^ 

- (y~ 1)(9 ^ . ^(f) ) 


T = 1 
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P = 1 — 

3^^<{>“V^<j) (1 - (Y~l)^^'i>) + 2 V(j) • V ( 3^ (1> ) 

+ (1 - (y- 2) a (j))w = 0 . (2.26) 

For a cylindrical combustor, it is desired to 
investigate the stability of the steady-state solution of 
(2.26d). To do this the steady-state solution must be 
obtained. It is assumed that the steady-state burning rate 
depends on z only (w = w(z)) and that the flow is axial 
(<j) = ±Cz)). Thus (2.26) simplifies to 

du/dz = w (]^ - d^/dz) 

T = 1 - hCy~l)u^ 

P = 1 - 

£ = 1 - ?gu2 . (2.27) 

It can be seen from (2.27) that the deviations of the 
steady-state solution from a uniform state are O(u^). To 
allow a discussion of orders of magnitude of various terms 
define the ordering parmeter e such that 


e = max (u) . 


( 2 . 28 ) 
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This implies that 

u = 0(e), w = 0(e), _^ = 0(e). (2.29) 

The stability analysis will now be carried out by assuming 
that 

<j> = K z) + e<}> ' ( X, y, z,t ) 

w = wCz) + e^w ' (x, yj z,t ) (2.30) 

where <fi' = 0(1), w’ = 0(1). (2.31) 

It is being assumed that the unsteady perturbation of the 
velocity potential from the steady state is of the same 
order of magnitude as the deviation of the steady state 
from a uniform state and that the unsteady burning rate 
perturbation is of the same order as the unsteady nonlinear 
gas-dynamic terms. The first assumption is necessary to 
be consistent with the quadratic approximation inherent in 
(2.26) while the second eliminates nonlinear terms Involving 
products of w' and (j> ' . These assumptions are made for 
simplicity and are not meant to imply that other orders of 
magnitude for the various terms are impossible or even 
unlikely. The objective of this work is to pick one case 
which is mathematically tractable and subject it to extensive 
analysis. While it is felt that most of the features of the 
stability problem will be reflected by this case, it is 
recognized that many other possibilities remain to be 
explored . 
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Substituting ( 2 . 30 ) and 

p = £ + ep ’ , P = P + eP’ , T = T + eT' (2.32) 

{p* = 0(1), P* = 0(1), T* = 0(1)} Into ( 2 . 26 ), retaining all 
terms of 0 (e'^) or lower, and dividing all resulting equations 
by £ leads to 

- V2({,’ + 2u9^^())’ + 

+ e { ( Y“1 ) + 2V4>**V(9^^’) +w’} = 0 

p' = _ { 9 ^ 4 ,. + u9^<),' + 3se(V())’ -V.}.' - (2-Y)(8^<i)’)^)} 

P' = _ y{9 4’ + U9 (|)’ + heCV4>' - (9^(fi’)2)} 

t — z t 

T’ = =Cy-1)(3^<{.‘ + u 9^^’ + 5seV<}.’ -V4.’). (2.33) 

The equations derived from this perturbation analysis 
are expected to be valid as long as the amplitudes of the 
perturbation quantities are finite and smaller than unity. 

For simplicity the primes appearing In the above 
equations will be dropped and It will be understood that 
unprlmed quantities are associated with perturbation from 
the steady state. 



Chapter 3 


ANALYTICAL SOLUTION OP WAVE EQUATION 

In this chapter a second-order solution for the 
velocity potential associated with transverse acoustic wave 
motion of an Ideal gas in a circular cylinder is obtained. 
The solution is found in the form of a Fourler-Bessel 
series. This provides a standard to compare the proposed 
methods which will be discussed in the following chapters. 

The governing equation is obtained from (2.33a) by 
assumption that the burning rate function is absent. Then 
the equation reduces to the second-order form of the one 
discussed by Maslen and Moore (2) for pure gas motion in a 
cylinder. 

9^^<j) - V^(j) + eC2V(j) • V 3^<j) + Cy-l) V^(j) 3^(j) ) = 0 (3.1) 

For a cylinder of radius. L^, a set of cylindrical 
polar coordinates (rj0jz) with the z axis being coincident 
with the cylinder's axis of symmetry is used. In seeking 
a second-order transverse wave solution, one assumes 

()) = (() (r,0,t) + E(|) (r,0,t) + ... . (3.2) 

1 2 

Substituting (3.2) into (3*1) and (2.33) leads to 
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= 0 

tt 1 1 

- v2(j) = _ (a,(u^) + (Y-I)v^(j) 3,(|) ) (3-3) 

^^2 2 1^1 

P = 1 + eP + e2p (3.4) 

1 2 

where Pj = - Y3,<t> , P = -y(9+-^ + - (9 4 , )2)) (3-5) 

tl2 ^2 1 

V2<j,^ = + 9^c^./r + 3Qe<)>i/r2, i = 1, 2 

u2 = (9 4) )2 + C9o<{> /r)2. (3«6) 

1 1 

Three different cases are discussed below. For the 
initial conditions 

Ki’jejO) = J (S r)cos0, 3, (|)(n,0,O) = 0 (3*7) 

111 ^ 

the solution of C 3 * 3 a) is 

4 > = J (S r)cos0cos(S t) 

1 1 11 11 

Substituting this equation to (3- 3b) produces 

9+-i-(!) -V^4> =(b + E E J^(S„v,r)cos (m9 ) )sin(2S t) (3-8) 

1 

where J is used to denote an n’th order Bessel function of 
n 

the first kind, is used to denote the m’th root of 

the equation J’ =0, and a prime indicates differentiation 
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with respect to the argument of the Bessel function. The 

coefficients b are integrals of Bessel functions which 
mn ^ 

must be computed numerically. The details of this 
calcuatlon are discussed in appendix A. 

Solving C3«8) one obtains 

<j) = (b /(4S^ ))(2S t - sln(2S t)) 

2 00 11 11 11 

■'Jo nil ‘ 'IS^ ^ )J„(S„nr)oos (me ) t5ln(2S ^ ) 

m7^ 1 

- ( 3 - 9 ) 

The other two solutions to be discussed subsequently are 
found in the same way. For the initial conditions 

(})(r3e,0) = Oj 3,<{>(r,0,O) = S J (S r)cos0 (3.10) 

t 11111 

it is found that 

(j) = J (S r)cos0sln(S t) (3.11) 

1 1 11 11 

and that is given by the negative of (3.9). For the 
initial conditions 

(})(rj0,O) = (S^^r )cos0 , 

3^(j)(r,0,o) = S J (S r)sine (3.12) 

^ 11 1 11 


the solution is 



2Q 


)cos (S 2 _ 2 _t -6) 

and <t> ^ t 2 b 2 j^/(S|^ -4s^^) }J 2 (S 2 j^r ) ( sin( 26 ) {cos (2S^^t ) 

- cos(S 2 nt)} + cos(20){sln(2S^3_t) - (2S^^/S2^)Sin(S2^t ) } ) . 

It can be seen that for these initial conditions 
represents a spinning wave but (j)^ + e ^2 does not. Some 
typical results for the pressure functions and P 2 were 
computed using the solutions for initial conditions (3.7) 
(labeled standing wave) and initial conditions (3- 12) 

(labeled traveling wave) for r = 1, 0 = 0, and y = 1.2. 

These data are presented graphically in Figure 1. P'or 
these computations the series expansions were terminated at 
n = 5- There is virtually no difference between the results 
for n = 4 and n = 5, and even the results for n = 1 provide 
a reasonably accurate solution. 

The results discussed above were used to check the 
numerical calculations to be discussed in the next two chap- 
ters. These results generalize those of Maslen and Moore (2) 
to Initial conditions which do not lead to periodic solutions. 



Chapter 4 


COLLOCATION 

In previous investigations of velocity-sensitive 
combustion instability {see, for instance, (6) and the 
references therein} the most widely used burning-rate 
equation is that associated with the vaporization limited 
combustion model discussed by Priern and Guentert (3). In 
the notation of the present thesis the second-order version 
of this equation can be expressed as 

w = Cw/e^)C(l - ^e8^(j))(Cl - u^/u^^) 2+(u^/u^) ^ )'^-l ) (4.1) 

where u is the component of perturbed gas velocity in the 
axial direction and u^ is the perturbed component normal to 
the axis of symmetry, and u^ is the magnitude of the droplet 
velocity vector which assumed to be axial. Substituting 
(4.1) into (2.33a) yields 

+ w9^(J) + 2u3^^(J) - V^(j)(l - e(Y~l)^^<i>) 

+ 2e^(j)"^9 <j) + Cw/e)((l - %e9, i|))((l - 3 (j>/u^)2 

t ~ t z n 

+ Cv<p-V<j> - C 9^<|) ) 2)/u2 )^- 1) = 0 . (4.2) 
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No exact solutions are to be expected to (4.2). Fully 
numerical solutions, on the other hand, are expected to be 
quite time consuming for multidimensional problems. For 
these reasons it was decided to attempt approximate 
partially analytical solutions using the method of weighted 
residuals {see for instance, (6)}. For problems of pressure 
sensitive combustion instability Powell (1) successfully 
employed the Galerkin methods This method is limited in 
applicability, however, to equations containing nonlinear- 
ities Involving only polynomial functions of the dependent 
variable and its derivatives. It is clear that (4.2) is not 
of this form. The orthogonal collocation method was, 
therefore, chosen because of its simplicity and adaptability 
to equations containing nonlinearities of any algebraic 
form. The application of this method to the problem of 
transverse wave motion in a cylindrical chamber will be 
discussed below. 

Consider the transverse wave motion in a cylindrical 
combustor. It is convenient to describe this problem in 
terms of cylindrical polar coordinates (r,0,z). Then 

<|) = (|)(r,e,t). (4.3) 

It will also be assumed that w and u^ are constants. Then 
(3.2) becomes 


9 <}>■*■ (() — ( 9 <f) + 9 (t)/r + 9 (j)/r^)Cl ~ e Cy~1 ) ) 

tt~t rrr ^ 
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+ 2e(a <f> 9 , (J) + 8 <()9 <|)/r2) + (w/e)((l - ^e9,(j))((l 

X* x*0 o dZ x> 

+ C(3^<{.)2 + (3Q4,/r)2)/u2)^ - 1) = 0 . (4.4) 

Now, assume an approximate solution of the form 




P Q 
z z 
m= on=i 


J (S r)cosCme )f (t) . 
m mn mn 


(4.5) 


Equation (4.5) is a superposition of the normal modes 
associated with the corresponding linear acoustic problem. 

A solution of this form allows one to Investigate the 
Influence of nonlinearities on the behavior of modal 
amplitudes which would exhibit simple harmonic motion in 
the linear case. Only standing modes can be described by 
(4.5). To represent spinning modes another series involving 
sin (me) must be added. This is omitted for simplicity in 
this discussion. It is convenient to rewrite (4.5) as 


<!> 


N 

Z J 


m 


J 


(S r)cos(m e )f (t) 
ni.ri. J J 


where N = (P + 1)Q . 


(4.6) 


The equation is now evaluated at N points (the collocation 
points) to yield 


4>i 


N 

Z C 


IJ^j 


1 


1,2 , . . .N 


(4.7) 
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where 


V 0 J 


(4.8) 


Next, (4.7) is Inverted to obtain 


N _i 

f. = S C. 7()). 

1 j=i iJ J 


(4.9) 


Then the various derivatives of (4.6) can be expressed in 

terms of value of d) . as 

1 


N 


N 


X ' 1 * 2 ? 

(ad)) = c .q, , = .z c (|) 

r 1 .i = 1 0 1 j = 1 ij J 


N 


N 


(a„<j)). = z C” <f, C9«n<t>). = 1=1, ...N (4.10) 


La , ill , 

j=l J 


'60^ 'i 


D = 1 


^ij "J 


where 


cj’. = ( 9 c. Jc cy^ = ( a c )C t 
ij kj " ij r.r. ik kj 

1 11 


-1 


c®. = (a.c.,)crl, c?® = (a „ c., )cr^. 

ij ®^ik kj " ij 0^6^ ik^ kj 


(4.11) 


Substituting into (4.4), a set of N ordinary second-order 
differential equations having the following form is 
generated. 


N 


N N 


+ W4>, - E C:,<|>,(l-e(Y-l)<f, ) + 2e E E C, ,.<(>,<}) 


j = l 




1:1 l,7i-ijk'^j'*'k 


N N 

+ (w/e ) ( (l-5ge(j) . ) (1 + E E C. , 4i , 4>. /u^ ) ‘*-l)=0 (4.12) 

- 1 j = i k=l ijk j k L 




where 
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(4.13) 


Solutions were obtained by a fourth-order Runge-Kutta 
method to the set of equations (4.12). It was found that 
the results were very sensitive to the number and locations 
of collocation points, especially to the radial distribution 
of points. It was found that even in cases when w was equated 
to zero (no combustion) it was possible to find many choices 
of collocation points which lead to the computed modal 
amplitudes increasing without bound. 

In order to illustrate the difficulties involved in 
a simple case, consider a one-dimensional problem of 
longitudinal wave motion with y = 1, u = 0 and w = 0. In 
this case (4.2) reduces to 


- 8 4 . + 2 eO 4>3 = 0 . (4.l4) 

tt zz z zt 


First consider the boundary condition 


3 4 .( 0 , t) = 0 , 4 >(T:,t) = 0 . (4.15) 

z 

A one-term solution satisfying (4.15) is 

4 . = f (t ) cos (^sz ) . (4. 16 ) 

Applying the collocation method with a collocation point zq 
produces an equation for 4 >g = 4 >(zQ,t) of the form 
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(}> + it;A + JsiCtan^Cz /2)d> d) ) = 0. (4.17) 

0 0 0 0 0 

The Galerkln method (to be discussed in detail in the next 
chapter) when applied to the same problem yields an equation 
for f of the form 

f + 4if + 4 ff/(37r) = 0. 

As another example consider the boundary conditions 
9^<t)C0,t) = 0, 32<f>CTr,t) = 0. 

A solution satisfying (4.19) is 

= f(t)cos(z) 

for which the collocation method leads to 

(!>+(()- 2etan2z <}> c|) =0 (4.21) 

0 0 0 0 0 

and the Galerkin method leads to 

f + f = 0. (4.22) 

Prom the above two cases, one can observe that the 
last term in each of the equations obtained by collocation 
can take on any value between 0 and “ depending on the 
location Zg of the collocation point. No such difficulty 
is encountered when using the Galerkin method. Even if more 
points were used, the same behavior was obtained. Therefore 


(4.18) 


(4.19) 


(4., 20) 
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there appears to be no way to rationally decide which sets 
of collocation points can produce the correct result. For 
this reason, after the expenditure of much effort, the 
collocation method was abandoned and the Galerkln method 
was adopted. This will be discussed In the next chapter. 



Chapter 5 


THE GALERKIN METHOD 

The Galerkin method is a special application of the 
method of weighted residuals (usually referred to as MWR) . 

It has been extensively used in the solution of various 
stability and aeroelastlcity problems {see for Instance 
Powell (1) } and proved Itself as a useful tool for the 
solution of both linear and nonlinear problems. Although 
it is an approximate mathematical technique, it has 
nevertheless produced results which were in excellent 
agreement with available exact solutions. These approximate 
solutions are usually simpler in form than the exact 
solutions obtained by numerical integration, and their 
quantitative evaluation requires considerably less 
computation time. However, this method is known to be 
reliable and applied conveniently only to equations 
involving nonlinearities of a polynomial type. Therefore, 
if this method is used in conjunction with the varporlzation 
limited burning rate function of (4.1) the problem becomes 
intractable. Thus, the following simpler purely phenomenolog 
ical burning-rate function was employed for the case of 
instantaneous combustion response. 

w = wnu^ (5.1) 


28 



29 


where n is a constant which will subsequently be referred to 
as the interaction coefficient. It can be determined either 
by comparison of predictions of the theory directly with 
experiment or by comparison with burning rate laws meant to 
apply to special types of combustion processes. As an 
example of the latter method consider the vaporization- 
limited burning rate law (4.1). As it stands (4.1) exhibits 
both pressure and velocity sensitivity. A purely velocity 
sensitive law can be obtained by assuming e << 1 to get 

w = (w/e^)((l + (u/u. )2)^ - 1) (5.2) 

V — i-i 

where the v is used to denote the vaporization model. If one 
equates, the slopes of (5.1) and (5.2) at u^ = 0 and plots 
both functions, the graph would look as follows. 
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It can be seen that (5.1) would overestimate the burning 
rate. Thus the upper bound for n can be computed. Prom 
(5.1) and (5.2) 

'^^2^ “ ^u2^V Cw/(4e2u2) )/(l+(u/u^) 2)3/^. 

Thus d ,w(0) = wn, d ,w (0) = w/(4e2u2). 

U'^ — U*^V — J_i 

Equating these results, one gets 


n 

max 


l/(4e2u2). 


C5.3) 


The phenomenological law C5.1) can be related to other 
special burning-rate laws in a similar manner. 

If it is desired to consider history-dependent 
combustion processes ^this can be done through the general 
formula 

t 

w = n/ G(t-?)d (u2)dg (5.4) 

0 5 

where G is memory function and C is a dummy variable. If 
G(t) = H(t) (H being the unit step function) all Increments 
of change in u2 occurring in the past are counted equally 
and (5.1) is recovered. If G(t) = H(t) - H(t-x) all 
increments of change in u2 are counted equally up to x units 
of time in the past while those previous to that time are 
not counted at all. Substituting this expression into (5.4) 
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and integrating one obtains 

w = wn(u2 - u^) (5*5) 

_ ^ 

where u =uCt-x) and t is called the time delay. Equation 

T 

(5.5) will be used to account for the history of the burning 
process in a rough way. More sophisticated treatments of 
this phenomenon could be obtained by employing a more 
realistic function for G(t). Substituting (5.5) and (5*1) 
into (2.33a) one obtains 

+ 2u 3^^(|) + wS^cj) + e( (y-l)V^(pd^(p + 2V(j>-Vd^(t> 

+ wn(V(})*V(j) - jV(f) «V(j) )) = 0 (5.6) 

— XT 

where j = 0 for instantaneous combustion and j = 1 for ; 
history-dependent combustion. 

The most general solution of (5.6) (subject to hard- 
wall boundary conditions for the unsteady variables) can be 
written in the form of the following Pourler-Bessel series. 

+ E E (f (t)cos(me)+g (t)sln(me))J (S r) (5*7) 
m=in=i mn mn ran 

Substituting (5-7) into (5-6) and applying the usual 
Galerkin orthogonalizatlon procedure leads to an infinite 
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set of coupled ordinary differential equations governing the 

functions f and g_„. No solution can be obtained unless 

mn mn 

the series (5-7) is truncated so as to produce a finite set 
of equations. If the terms neglected are actually small, an 
accurate solution will result. 

In this thesis attention will be focused on initial 
disturbances having the form of the first tangential mode. 
The simplest finite series contained in (5.7) capable of 
modeling the effect of quadradic nonlinearities in (5.6) is 

= f^(t )J^(S^^r) + f^(t )J^(S^^r)cos0 + f ^(t ) J^(S^^r )cos26 

+ g2(t )J2_(S^^r)sin0 + g2(t )J2(S22^)sin20 • (5.8) 

The Galerkln method then produces the following ordinary 
differential equations governing these five variables. 




+ CjCfgfj + gjg^) + + gp + C^3(f| + gp 


-J('=lTox " ° 


q + + CgCg^g^ f fpp 


+ cpg^g^ + fpp + wn(C3^fp^ + + g^gp 
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^2 “'' -^2 '*' ^21^2 ^9^2^0 '*' ^10 


+wn(C^g(f2 g^)+C^^fQf2~ S^^) + ^i7^0t^2t^^^ ^ ° 

§1 + ES^ + ^ii^i ■'■ "'' '^6^^1®2 “ ®1^2^ 

+ C^(g2i\ - f2gi) + wnCC^^f^g^ + “ ®1^2^ 




^2 ^^-^2 ^ 21®2 ^ ^^" 8 '‘ 0®2 S ^ 2"'0 ~ "lO Vl 


+ f^S^) + wnCC^^f^g^ + 2C^gf^g^ 




(5.9) 


The coefficients C^ (1=1,2, ... 17) are integrals of 
Bessel functions which must be computed numerically. The 
details of this calculation are discussed In appendix A. 

The stability boundaries in the (n,e) plane 
were determined first. For a given value of w, a value of 
e was selected and solutions were obtained for various values 
of n. In each case, if the modal amplitudes exhibited 
growth with time, the value of n was decreased for the next 
run while if decay of modal amplitudes was observed, the 
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value of n was Increased accordingly. A systematic 
iteration process was devised so that the computer could 
carry out these calculations without analyst Intervention. 

In this wayjOne point on the stability boundary was 
established. Then a new value of e was chosen and the 
iteration process was repeated to determine another point 
on the stability boundary. This procedure (which consumed 
a considerable amount of computer time) was continued until 
enough points had been found to establish the shape of the 
stability boundary. It should be pointed out that the 
amplitude of the IT mode always initially decreases due to 
the fact that purely velocity-sensitive combustion 
instability is always linearly stable. Thus it is necessary 
to continue the calculation for a considerable period of 
time to determine whether a given set of conditions 
corresponds to nonlinear stability or instability. 

Figures 2 and 3 show some typical stability 
boundaries for instantaneous combustion using the initial 
conditions 

f^(O) = 1, f^CO) = fgCO) = g^(0) = ggCO) = 0 
fo(O) = ^1^0) = ^2^0) = Sl^°) = S2(0) = 0 (5.10) 

and f^(0) = 1, f^CO) = f^CO) = g^(0) = g^CO) = 0 

giCo) = 1, f^co) = f^(0) = f^co) = g^co) = 0, (5.11) 
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respectively. The region below and to the left of the 
stability boundary is associated with stable conditions, 
while the region above and to the right is associated with 
unstable conditions. 

In the case of linear acoustics, the first set of 
initial conditions would lead to a standing wave and the 
second set of initial conditions corresponds to a traveling 
wave. In both cases, it can be seen that the stability 
boundaries have roughly the forms of rectangular hyperbolas. 
As expected, increasing the steady-state burning rate reduces 
the value of n required to produce instability. This effect 
is somewhat more pronounced for the traveling waves than for 
the standing waves. There does not appear to be a distinct 
pattern to the results. Thus for w = 0.2 standing waves are 
more stable than traveling waves while for w = 0.1 traveling 
waves are more stable than standing waves. 

Figures 4 and 5 show stability boundaries associated 
with initial conditions (5.10) and (5.11) > respectively, for 
history-dependent combustion. It is apparent that the 
Influence of the time delay parameter on the results for 
standing waves is much greater than on the results for 
traveling waves. Again, no clear pattern emerges from the 
results. Comparing Figures 2 and 4, and 3 and 5 shows that 
instantaneous combustion can be either more or less stable 
than history dependent combustion depending on the value of 
the time-delay parameter. Comparing Figures 4 and 5, it 
can be seen that standing waves can be either more or less 
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stable than standing waves for history dependent combustion. 
Figures 4 and 5 also Illustrate the fact that increasing the 
amount of time delay can either increase or decrease the 
stability of the system. Thus^for both standing and 
traveling waves, x = ir corresponds to a more stable situation 
than either x = 0.5it or x = 1.5ir . The lack of patterns 
exhibited by these results emphasizes the need for numerical 
methods of the type developed during the present research. 

It appears that the only way to find out what will happen in 
a given situation is to solve the equations for that 
particular case. This matter will be discussed further 
subsequently . 

The pressure perturbation can be calculated from the 
velocity potential perturbation by substituting equation 
(5.8) into (2.33c), expanding the result in a Fourler-Bessel 
series and retaining only the terms corresponding to the IT, 
2T and IR modes to obtain 




+ *^ 05^0 ^ Sl^ ^ 07^^2 ^))) Jo ^^ 01 ^^ 


((diifi + e(d^2^'o^l ^13^^1^2 


+ '^ 14 ^ 0^1 
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+ d^^(f^g 2 - f 2 S 3 _) ) )sln 0 )J^(S^^r ) 


+ ({d^ifj + =('^ 2 / 0 ^ + d23'fi - Si> + '524‘'of2 


• • 


+ d (f 2 - g2 ) ) )cos2e 
1 1 


+ (d2^g^ + e(d22fQg2 + 2d22f']_S2 + ^24^0^2 


+ 2 d 2 ^f^g^) )sln 20 )J 2 (S 2 ^r ) ) 


(5.12) 


The coefficients d„ , d^ . d^ are calculated in the same 

Os Is 2s 

fashion as those discussed previously and the values are 
listed in table 8 of the appendix. 

Some typical results for wall pressure waveforms are 
presented in figures 6 through 37* 

Figures 6 through 9 correspond to instantaneous 
combustion with e= .05, w= .1, n= 175, and initial 
conditions (5.10). This leads to a stable standing wave 
oscillation and the pressure can be observed to decrease 
gradually. However ,by changing the interaction index to 
n = 220 (Figures 10-13) the pressure is seen to grow 
gradually. This is an unstable situation. In both cases 
the response is dominated by the IT mode but distorted to 
some extent by the presence of the 2T and IR components. 

Figures 14 through 17 correspond to history -dependent 
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combustion with n = 3 ^ 5 , while the other parameters are the 
same as before. This is a stable situation and the amplitude 
of the pressure slowly decays with time. 

Figures 18 through 21 are obtained by changing the 
interaction index to n = 352. This is an unstable situation 
as indicated by the growth of the pressure amplitude with 
time. In both of these cases the presence of the 2T and IR 
components is much more noticable than it was in the first 
two situations. 

Figures 22 through 25 correspond to instantaneous 
combustion with n = 180, e = .05, w = .1, and initial 
condition (5-11). This produces a traveling wave. The 
pressure is observed to be decreasing with time (a stable 
situation) while changing the interaction index n to 210 
(Figures 26 through 29) causes the pressure to increase 
(an unstable situation). Figures 30 through 37 show the 
significance of the time delay function. Figures 30 through 
33 are obtained by setting t = tt and n = 197. This is a 
stable case. The pressure is observed to decrease. Setting 
n = 199 (Figures 34 through 37), on the other hand, produces 
an unstable situation when the pressure increases very 
rapidly. In both of these situations the response appears 
to be dominated by 2T contribution to the pressure. 

From these figures one can conclude that the pressure 
waveforms exhibit a strong second harmonic distortion and 
this distortion arises from the effect of the quadratic 
nonlinear terms. A variety of behaviors are possible 
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depending on the nature of the combustion process and the 
parametric values involved. 



Chapter 6 


ONE-DIMENSIONAL MODEL 

In this chapter an annular combustion chamber with a 
gap width much smaller than the Inner radius Is considered. 
This geometry will henceforth be referred to as that of a 
narrow annulus. While this geometry Is not of much practical 
Interest, It Is quite useful for the purposes of analysis. 

This Is because only one space coordinate Is needed to 
describe the problem.- This makes a direct finite-difference 
numerical solution of the original partial differential 
equation feasible and also simplifies the algebra required, 
to carry out the Galerkln modal analysis. In what follows 
three questions will he Investigated. First, the effect of 
changing the numerical values of certain coefficients 
appearing In the governing equations for the modal amplitudes 
will be discussed. Second, the Galerkln solution will be 
checked using a finite difference numerical solution of the 
complete equation. Third, numerical solutions of the 
complete wave equation using the vaporlzat lon-llmlted 
burning-rate law will be compared to similar solutions 
associated with the phenomenological burning-rate law 
employed In the previous chapter. 

The appropriate wave equation for transverse wave 
motion in a narrow annulus can be obtained from (5.6) by 
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setting r = 1 and 3^ = 0 . To further simplify the results 
the parametlc values y = 1 (isothermal process) and J = 0 
(instantenous burning response) will be employed. For this 
situation (5.6) simplifies to 

+ 2e3g(|)3g^(j> + wn e(3g(j))^ = 0. (6.1) 

To carry out the modal analysis, it is assumed that 
the potential function can be expressed as 

(j) = f^(t)cos0 + f2(t)cos29 + g^(t)sine +g 2 (t)sin 2 e. (6.2) 

Then applying the usual Galerkin orthogonallzation 
procedure leads to 


fl + G^f^ + wf^ + eA^(f^f 2 + ^ 2^1 S3_g2 + § 261 ' 


+ ewnA (ff +gg)=0 
- 2 1 2 12 


^2 ^2^2 -^2 ''' + £wnAi^(g| - f|) = 0 


Si + afSi + wEi + + fiSa - e/s - VF 


+ emik^ir^S2 - ^ ° 


&2 ^ 2^2 -^2 ” ^^3^®1^1 ^l^i) “ 2ewnAi|f^g2^=0 (6.3) 
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where = 1, = 2^ = 2 , - '2- ^ k^ = 1 and = .5. 

The symbols A^, A 2 , A^ and k^ have been inserted to 

illustrate the effect of changing their numerical values. 

McDonald (7) in his analysis of combustion- 
instability in an annulus found that for a given value of e 
the value of n required to produce instability was approx- 
imately twice as high for standing waves as for traveling 
waves. The results discussed in the previous chapter for 
a full cylinder indicated no such relationship. Equations 
( 6 . 3 ) were employed in an attempt to determine the factors 
which have a significant effect on the stability boundary. 
Several calculations were made and a representative sample 
of the data thus obtained is presented in Tables 1 and 2. 

It was determined by McDonald that the terms representing 
gas -dynamic nonlinearities, had a small qualitative effect 
on stability calculations. Thus the coefficients A^ and 
vjere held fixed during these calculations. 

The entries in the first two lines of each table 
were computed using the correct equation for an annulus. It 
can be seen that the standing wave is twice as stable as the 
traveling wave, in agreement with the results of McDonald. 
The third and fourth lines in each table were computed by 
changing A^ from .5 to .77 (This makes the ratio A^/A^ 
the same as the ratio of the corresponding terms in the 
governing equations for the full cylinder.). It can be 
seen that this lowers the stability limit in all cases but 
does not alter the fact that an initial disturbance in the 
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Tab le 1 


Stability Limits for Standing Wave 
in Annular Combustor 


n 

^2 

A 4 

e 

91 

2 

.5 

.05 

46 

2 

.5 

.1 

74 

2 

.77 

.05 

37 

2 

.77 

.1 

735 

1.66 

.5 

.05 

379 

1 . 66 

.5 

.1 

595 

1.66 

.77 

.05 

30 7 

1.66 

j 

.77 

.1 


Table 2 


Stability Limits for Traveling Wave 
in Annular Combustor 


n 

S^2 

A 4 

e 

44 

2 

.5 

.05 

22 

2 

.5 

.1 

36 

2 

.77 

.05 

18 

2 

.77 

.1 

66 

1. 66 

.5 

.05 

31 

1.66 

.5 

.1 

54 

1.66 

.77 

.05 

26 

1.66 

.77 

.1 
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form of a standing wave is twice as stable as one in the 
form of a traveling wave. 

The fifth and sixth lines in Tables 1 and 2 are found 
by restoring to its original value .5 and changing 
from 2.00 to 1.66 (This makes the ratio $^ 2 / 0 ^ equal to the 
ratio 322 / 32 ^ 2 . 3-ssociated with the full cylinder,). For this 
situation it can be seen that the standing wave is approx- 
imately ten times as Stable as the traveling wave. Further- 
more , the effect of the value of JJp 4he location of 
stability boundary is much greater for a standing-wave 
motion than for a traveling-wave motion. The last two lines 
are associated with implementing both of the changes dis- 
cussed above simultaneously. These results confirm that the 
change in lowers the stability limit in all cases but 
does not affect the relative stability of standing -wave and 
traveling-wave disturbances. 

The data presented above indicate that standing waves 
will be twice as stable as traveling waves only under very 
special circumstances (^ 2 /^ 2 . “ There is no reason to 

expect this to be a characteristic of other systems exhibit- 
ing combustion instability as, in fact, it is not for a full 
cylinder . 

In this section, a finite different procedure is 
employed to generate a set of second order differential 
equations. The central difference formulas 

'a'*’! = 


)/C2Ae) 

x+l 1-1 





(6.4) 


are used to produce the following set of governing equations 


(j)^ = -w<()^ + “ 2(})^ + (j)^_^)(l - e(Y-l)i}>j_)/(Ae)^ 

- - <1>._3_)/(2A0)2 

- ew^ 2<i<N-l (6.5) 


where N is the number of points employed. A forward 
difference formula 


^ e'*’l " ^“^^1 ^ ^'**2 -^ 3 )/ C 2 A 0 ) ( 6 . 6 ) 

and backward difference formula 

^0'^N ~ 2 ~ ^‘^N— 1 (6.7) 

are used for the boundary conditions 

3q<1)(0) = 0, 9q()>(tt) = 0. (6.8) 


Hence, along the boundaries, (6.5) becomes 



For the phenomenological model,. the turning rate 
function becomes 


w^ = nw(<t,^_^^ - (}.^_^)2/C4 a 92) 2 < i < N-1. (6.10) 

Combining (6.8) and (6.10) one obtains 

W2 = 4nw(4>3 - (f>2) ^/(9A0 ^ ) 

^N-1 " ~ (6.11) 


respectively . 

For the vaporization model, the burning rate 
function becomes 


= w((l -e^^/2)(l + 

- (j). T )2/(2A0n )2)^ _ i)/e2 2 < i < N-1 (6.12) 

X"" X X/ 


and 3 along the boundary, yields 
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w - w((l - ) (1 + (2(4) - 4>o ^/^3A0 u ) ) ^ ) “* - l)/e^ 

d. j d L 

accordingly . 

Using either the phenomenological or.,, the yaporlzatlon- 
llmited combustion model, (6. 5) and (6.9) can be solved by 
the Runge-Kutta method to obtain cj)^, 1 = 2,3j...j N-1. 

Then the modal amplitudes can be obtained by using the 
Pourler-cosine transform 


IT 

f.(t) =2/ (},cosi0de/TT. (6.14) 

1 0 

These Integrals must be computed numerically since (j> is 
known only at the grid points. 

Some typical results for standing waves are presented 
In Table 3- They are Intended to Illustrate the 
accuracy of the two-term Galerkin solution and to compare 
the results associated with the phenomenological and 
vaporization-limited combustion laws. The column labeled 
PG indicates the results with a two-term Galerkin solution 
using the phenomenological model, the column labeled PP 
contains the results of a finite-difference solution of 
these equations, and the column labeled VP presents results 
of a finite difference solution using the vaporization- 


limited model. 
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Table 3 


Stability Limits For Phenomenological and 
Vaporization-Limited Combustion Model 



PG 

PF 

VF 

e 

n 

n 

^L 

.1 

46 

28 

1.3 

.05 

91 

56 

2.1 

.01 

452 

286 

5.2 


It can be seen that the order of magnitude of the 
interaction index n required for Instability for both the 
finite difference and Galerkin methods are roughly the same 
but the stability boundaries predicted by the Galerkin method 
are approximately twice as high as those predicted by the 
Finite-difference method. It is to be expected that the 
Galerkin method will lead to higher stability limits than 
the use of an exact solution procedure. This can be explained 
as follows. The instability mechanism is basically a feed- 
back process. Consider the form of (6.3) associated with 
standing waves in an annulus. This is 

fl + wf^ ^1 2eCfj^f*2 + ^ 2ewnf^f2 = 0 

^2 —^2 ~ ^^1^1 ~ %ewnf^ = 


0 . 


(6.14) 
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For clarity -the gas-dynamic nonlinearities will be neglected 
to give 

« 

+ wf^ ^1 2eTmf^f2 = 0 
+ wf^ + 4f^ - isewnf^ =0. (6.15) 

A one-terra Galerkin solution would lead to the 
equation 

f2_ ~ 

This would indicate unconditional stability. Instability 
can arise only when the presence of the last term in (6.15b) 
causes to grow and this then causes f^ to grow because 
of the last term in (6.15a). The growth of also produces 
the growth of higher modes which in turn causes additional 
growth of f^. The contributions of these higher modes are 
neglected in a two-term Galerkin analysis and thus the energy 
input due to unsteady burning is underestimated. For this 
reason the Galerkin analysis will overestimate the stability 
of the system. This overestimation will decrease as the 
number of terms retained is increased. An example of this 
can be seen by comparing the one- and two-term analyses . 

A one-term solution predicts that the system is always 
stable. A two-term solution predicts the correct qualita- 
tive behavior but overestimates the system’s stability by 
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rougb.ly a factor of tvfo. It is to be expected that further 
Increases . in accuracy can be obtained by keeping more terms 
in the assumed solution but that this will not seriously 
affect the quatltatlve prediction. 

In Chapter 5, it was found that equating dw/du^ the 
phenomenological and vaporization-limited combustion laws 
at u^ = 0 lead to the formula 

n = l/(4e^u^). (6.16) 

Assuming that the actual data can be fitted to a formula of 
the form 

n = C/(4e2u2) (6.17) 

the data give in Table 3 provide an opportunity to estimate 
C. The results are shown below* 


e 

.1 

.05 

.01 

C 

1.89 

2.47 

3.09 


It appears that a value of C = 2.5 would give accept- 
able accuracy over this range of e. 

For a given value of u^^ the n computed from (6.l6) 
will overestimate the burning rate. Thus it might be 
thought that C should be less than unity. The stability 
criterion was, however, based on the magnitude of the modal 
amplitudes. Inspection of Tables 4 and 5 shows that the 
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vaporization-limited combustion model Involves the higher 
modes to a greater extent than does the phenomenological 
combustion model. Thus, a given value of u^ can be associated 
with lower Individual modal amplitudes in the former case 
than in the latter. These two effects Interact and calcula- 
tion shows that C is greater than unity in this range of e. 
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TABLE A 

MODAL AMPLITUDE FOR PHENOMENLDGICAL MODEL 
IN ANNULAR COMBUSTOR 


T 

MODE 1 

MODE 2 

MODE 3 

MODE 4 

MODE 5 

MODE 6 

MODE 7 

MODE 3 

Q-CO 

1.000 

0.000 

-0,000 

0,000 

-0.000 

0 .000 

-0. 000 

0. 00 0 

1.25 

0.3 37 

0. 029 

0. 002 

0,000 

0.000 

-0,000 

-0.000 

-0.000 

2.50 

-0.687 

-0.013 

-0.000 

0.000 

o.oco 

-0.000 

0.000 

0,000 

3. 75 

-0.700 

0,0 64 

-0.003 

0.001 

-0.000 

0. 000 

0. 000 

-0.000 

5. CO 

0. 2C9 

-0. 043 

0.009 

-0.001 

“0.000 

0 .000 

-0.000 

0.000 

6.25 

0.751 

0.0 40 

-0.006 

-0.003 

”0.001 

0.000 

0.000 

0-000 

7.50 

0.237 

0.025 

0.000 

-0.002 

-0.001 

-0.000 

-0.000 

-0,000 

3,75 

-0. 535 

-0.060 

0.015 

0.004 

-0.002 

-0.000 

0.000 

0.000 

10. CO 

-0.523 

0.106 

-0.023 

0.009 

-0.003 

0.001 

-0. 000 

0.000 

11.25 

0.187 

-0.088 

0.024 

-0.000 

-0.003 

0.002 

-0.000 

-0.000 

12,50 

0. 585 

0.062 

-0.009 

-0.011 

-0.005 

-0.001 

0.001 

0. 001 

13.75 

0.153 

0,011 

-0.010 

-0, 010 

-0.006 

-0.003 

-0.002 

-0.001 

15. CO 

-0, AAO 

-0.065 

0.032 

0-007 

-0.007 

“0.001 

0.002 

0.000 

1 6.25 

“0.390 

0 .112 

-0.043 

0.020 

-O.OlO 

0,005 

-0. 002 

0.001 

1 7, 50 

0.183 

-0.105 

0.030 

0.005 

-0.010 

0 .004 

o.ooi 

-0.002 

1 8.75 

0.463 

0.078 

-0.004 

-0.017 

-0.013 

-0.004 

0.001 

0. 002 

20. CO 

0.085 

-0. 009 

-0.023 

-0.020 

-0.014 

-0,008 

-0. 004 

-0.002 

21, 25 

-0. 37A 

-0.050 

0.045 

0.004 

-0.015 

0 .00 1 

0.006 

-0.000 

22.50 

-0,287 

0 . 101 

-0.048 

0.028 

-0.013 

0,011 

-0.007 

0.003 

23.75 

0. 1 84 

-0.107 

0-025 

0.015 

-0.018 

0-005 

0. 004 

-0.005 

25. CO 

0. 371 

0 .089 

0.009 

-0,016 

-0.019 

-0,010 

0,001 

0. 00 5 

26.25 

0,030 

-0. C31 

-0.0,37 

-0, 030 

-0.020 

-0.012 

-0.006 

-0. 002 

27.50 

-0, 327 

-0.027 

0.052 

-0.005 

-0.021 

0,005 

0.010 

-0.002 

28.75 

-0.206 

0.081 

-0,044 

0.0 30 

-0,023 

0, 017 

-0, Oil 

0.006 

3 0. CO 

0, 1 91 

-0 . 102 

0,013 

0.028 

-0.Q23 

0.00 3 

0.008 

*“0 « 0 0 6 

31,25 

0.302 

0.093 

0,026 

-0.009 

-0,023 

-0.017 

-0.003 

0, 004 

32,50 

-0.019 

-0.054 

-0.051 

-0. 038 

-0,024 

-0.012 

-0-004 

0,001 

3 3.75 

-0.297 

0. COO 

0.056 

-0.019 

-0.024 

0.012 

0.010 

-0. 005 

35.00 

-0.141 

0.053 

-0.034 

0.027 

-0.024 

0- 0 20 

-0.013 

0.008 

36. 25 

0.207 

-0.095 

-0,005 

0.040 

-0.023 

-0.003 

0.012 

-0.005 

37.50 

0. 252 

0.108 

0.046 

0.004 

-0.023 

-0.022 

-0.008 

0. 001 

38,75 

-0,067 

-0. 079 

- 0. 0 64 

-0.042 

-0.022 

“0 .007 

0.002 

0,005 

AO. 00 

-0.284 

0.031 

0.056 

-0.036 

-0,020 

0.019 

0.005 

-0.008 

A1.25 

-C. 085 

0.033 

-0.019 

0.018 

-0.019 

0.016 

-0. 010 

0-007 

A 2. 50 

0,236 

-0.087 

-0.029 

0.052 

-0.017 

-0.013 

0.013 

-0.00 1 

A3. 75 

0,218 

0.120 

0.070 

0.023 

-0.014 

-0.022 

-0.014 

-0. 006 

A5. CO 

-0.122 

-0.109 

-0.076 

“0.039 

-0.011 

0, 006 

0-012 

0.010 

A6, 25 

-0.292 

0.0 70 

0.049 

-0.055 

-0.007 

0,025 

-0.007 

-0.008 

A7.50 

-0.029 

-0.001 

0.005 

-0.001 

-0.003 

0.002 

0. 000 

-0.000 

A3. 75 

0. 290 

-0.073 

-0.063 

0.058 

0.003 

-0.026 

0.007 

0.009 

50, CO 

0. 192 

0.135 

0.098 

0 .053 

0.011 

-0 .008 

-0.014 

-0. 01 3 

51,25 

-0,207 

-0, 152 

- 0. 0 81 

-0. 020 

0.016 

0.026 

0.019 

0.006 

52,50 

-0.328 

0.134 

0.022 

-0.073 

0.026 

0,015 

-0.022 

0.006 

5 3, 75 

0.055 

-0 .060 

0.055 

-0.044 

0.034 

- 0, 028 

0. 022 

-0.015 

55, CO 

0. 4 04 

-0.033 

-0.119 

0.037 

0.046 

-0.022 

-0.020 

0.012 

56, 25 

0,151 

0. 140 

0.125 

0.094 

0. 061 

0.038 

0.019 

0. 008 

57,50 

-0.403 

-0.219 

- 0. 0 45 

0.054 

0.065 

0.025 

-0.013 

-0. 022 

58.75 

-0.412 

0. 276 

-0.0 92 

-0.036 

0.071 

-0,051 

0.013 

0.012 

60.00 

0.324 

-0.246 

0.183 

-0.117 

0.060 

-0.017 

-O.Oll 

0,022 
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TABLE 5 

MODAL AMPLITUDE FOR VAPORIZATION MODEL 
IN ANNULAR COMBUSTOR 


T 

MODE 1 MODE 2 

MODE 3 

MODE 4 

MODE 5 

MODE 6 

MODE 7 

MODE 8 

0.00 

1.000 

0.000 

-0.000 

0.000 

-0.000 

0 ,000 

-0. 000 

0. 00 0 

1.25 

0.33A -0,193 

0, 0 39 

-0.007 

0.002 

-0.00 1 

0.000 

-0.000 

2-50 

-0. A55 

0,095 

-0.033 

0.009 

-0,002 

0,000 

0.000 

-0.000 

3,75 

-0.51A -0.111 

0.0 42 

0.C18 

-0.007 

-0,00 3 

0, 002 

0. 000 

5. CO 

0. 06 0 -0.035 

-0.036 

-0 .029 

-0 .020 

-0.012 

“0.007 

-0.003 

6.25 

0. AA3 

0.106 

-0.022 

-0.041 

-0.021 

-0,002 

0.006 

0.006 

7,50 

0.293 -0- 225 

0. 077 

0.020 

-0.043 

0.01 9 

0.007 

-0.013 

3.75 

-0. 32A 

0.160 

-0.104 

0.071 

-0.049 

0,034 

-0.024 

0.017 

1 0. CO 

-0.717 -0,029 

0.192 

-0.037 

-0.072 

0.042 

0,020 

-0.032 

1 1. 25 

-0.22A -0.A55 

-0. 330 

-0.172 

-0. 040 

0,041 

0,065 

0,047 

1 2.50 

1.281 

0.558 

0 o 3 0 0 

0.131 

0.040 

-0.003 

-0.021 

-0. 030 

13.75 

2.176 -0.505 

-0.474 

0. 08 7 

0,172 

0.024 

-0,078 

-0.05 3 

15-00 

0 .26A -0,967 

0.6 44 

-0.327 

0,068 

0.075 

-0.105 

0.057 

16.25 

-2.907 

0.8 57 

-0.311 

0.100 

-0.053 

0,005 

0. 019 

-0.031 

17.50 

-3. A16 -1.125 

0. 382 

0.037 

-0,193 

0.052 

0,072 

-0,074 

18,75 

0.AA8 -0.79A 

-0.551 

-0.324 

-0,057 

0.085 

0,110 

0. 056 

20.00 

3.619 

0. 5A6 

0. 200 

0.077 

0-001 

-0.024 

-0.039 

-0.044 

21, 25 

2.7A7 -1.259 

-0.224 

0.106 

0 • 166 

-0-024 

-0.0B7 

-0,025 

2 2.50 

-1.AA2 -0.328 

0- 340 

-0,271 

0. 145 

-0.042 

-0.032 

0.06 3 

2 3. 7 5 

- 3. 22 9 

0.136 

-0.113 

-0.068 

0,063 

-0,080 

0,058 

-0-038 

2 5.00 

-1.329 “1,0 66 

0.005 

0.181 

-0.028 

-0.102 

0 .012 

0, 059 
-0.049 

26.25 

1. 79A 

0, 024 

-0. 009 

-0. 085 

”0.086 

-0.083 

-0.067 

27,50 

2.388 -0.063 

-0.093 

-0,167 

-0, 089 

-0.030 

0.022 

O.OA 1 

28. 75 

0®^63 "0®881 

0,235 

0.090 

-0.114 

- 0, 006 

0. 069 

-0.027 

3 0. GO 

-1. 8C7 

0.272 

“0.195 

0.095 

-0.067 

0.045 

"0.031 

0.021 

31.25 

-2-038 -0.29A 

0,261 

-0.143 

-0.024 

0.075 

-0.052 

“0- GO 3 

32, 50 

-0.063 -0,796 

-0.404 

-0,100 

0,077 

0.093 

0.025 

-0. 040 

33.75 

2.115 

0, A17 

0.250 

0.120 

0.079 

0.049 

0.033 

0.021 

35.00 

2.172 -0.6A7 

-0, 302 

-0.013 

0-128 

0, 062 

-0. 032 

-0,059 

36- 25 

-0.AA7 -0.603 

0,422 

-0.241 

0.058 

0.049 

-0.079 

0,050 

37.50 

-2.586 

0. A3A 

-0.193 

0.045 

-0.004 

-0.025 

0.035 

-0. 039 

3 8.75 

-2.0AA -i.QOO 

0,179 

0. 126 

=0.127 

-0,044 

0.074 

0.000 

A 0.00 

1.1A2 -0.253 

-0.260 

-0.241 

-0.154 

-0.072 

-0.004 

0.038 

A1.25 

2. 3A8 

0.253 

0.055 

-0 .083 

-0 .089 

-o.oei 

-0. 054 

-0.027 

A 2. 50 

1. A97 -1,1A0 

0,032 

0.196 

0 .007 

-0.109 

0.006 

0.063 

A 3-75 

-1,721 

0.082 

0,011 

-0.091 

0.095 

-0.09 3 

0.077 

-0, 057 

A5. GO 

-2.759 -O.OIA 

0.116 

-0.188 

0.105 

-0.030 

-0.020 

0.044 

A 6S 25 

-0,807 -1,084 

-0,251 

0.126 

0.118 

-0.033 

-0.081 

-0.012 

A7.50 

2,050 

0, 301 

0.193 

0.083 

0,050 

0.024 

0,010 

-0.001 

A 8,75 

2. 51 A -0,289 

-0,260 

-0-186 

-0.011 

0. 062 

0.068 

0.023 

50,00 

0.206 -0.932 

0,405 

-0.043 

-0.117 

0,082 

0.014 

-0, 063 

5 1.25 

-2.255 

0. 397 

-0. 263 

0. 134 

-0.096 

0.066 

“0.051 

0.040 

52-50 

-2.287 -0.575 

0. 312 

-0.076 

-0.1C6 

0.084 

-0, 002 

-0.053 

5 3. 75 

0,331 -0.709 

-0,439 

-0.200 

-0,000 

0. 08 3 

0. 074 

0, 01 7 

55.00 

2.A77 

0.407 

0.229 

0.082 

0.033 

0.004 

-0.011 

-0.021 

56,25 

2.025 -0,858 

-0,236 

0.073 

0,141 

-0.002 

-0.072 

-0.031 

57,50 

-0,915 -0.377 

0. 322 

-0.252 

0.129 

-0.030 

-0.036 

0.060 

58.75 

-2.669 

0 . 303 

-0.109 

-0.033 

0.057 

“0.068 

0.057 

-0.042 

50.00 

-1.588 -1.060 

0,038 

0.179 

-0.051 

-0, 099 

0.028 

0.055 
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CONCLUSIONS 

The primary objective of this study was the 
development of a new analytical technique to be used in the 
solution of nonlinear velocity-sensitive combustion 
instability problems. Such a method should be relatively 
easy to apply and should require relatively little 
computation time. 

In an attempt to achieve this aim, the orthogonal 
collocation method was Investigated first. However, it was 
found that the results were heavily dependent on the location 
of the collocation points and characteristics of the 
equations. Therefore, the method was rejected as unreliable. 

Next, the Galerkin method, which has proved to be 
very successful in analysis of the pressure sensitive 
combustion Instability, was considered. This method proved 
to work very well. It was found that the pressure waveforms 
exhibit a strong second harmonic distortion and a variety 
of behaviors are possible depending on the nature of the 
combustion process and the parametric values involved. 

Finally, a one-dimensional model provided further 
insight into the problem by allowing a comparison of 
Galerkin solutions with more exact finite-difference 
computations . 


54 



55 


Some major conclusions of this research are as follows. 
(1) The form of (2.33a-) is somewhat insensitive to the 
specific assumptions used to derive it. This is demonstrated 
by the fact that Powell (1) developed an equation of a very 
similar form using a different set of assumptions. (2) The 
orthogonal collocation method is unsulted to solution of 
problems of the type under discussion here. (3) The Galerkln 
method is well suited to the solution of such problems. (4) 
Stability boundaries and pressure wave forms appear to be 
highly dependent on the parameters of the problem (interac- 
tion index, time delay, steady-state burning rate, etc.). 
Furthermore , there appears to be no clear pattern to the 
computed results. (5) The phenomenological’ burning-rate law 
employed in the majority of the work discussed predicts 
results which are qualitatively similar to those associated 
with the vaporization-limited burning-rate law used by 
previous investigators. The phenomenological law, further- 
more, can be used in conjunction with the Galerkln method 
while the vaporization-limited law cannot. (6) The computer 
programs developed in the course of this work can be 
employed to determine stability boundaries and pressure 
waveforms without the expenditure of excessive computer 
time. This is Important because the lack of clear trends 
discussed above make an individual analysis of each 


situation desirable. 
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Figure 1 


STANDINS WAVE 

Pressure functions for standing and traveling 
waves 














































































































APPENDIX A 


DERIVATION OP THE COEFFICIENTS 

In this section, the coefficients b... C and d . . 
appearing in chapter 3 and 5 are derived and their numerical 

values are presented in Table 6, J, and 8 respectively. 

In chapter 3, for initial condition (3-7), one assumes 

oo 

94-+.<J)o - = 2 H Cr)cos (me ) slnC2S- - t ) . 

w d. in= 0 ax 

By comparing coefficient yields 

HgCr) = inS^^CPSf^Jg - Cy-DSf^Jf - 4s^^jQj^/r + 4Jg/r2) 
H^Cr) = 0 

H2(r) = 3^S^^(2S2^J2 - (y-DsJ^jf - llS^^Igl^/r) 

H^Cr) =0 for m = 3, 4, . . . , 

Then expanding an a Fourier-Bessel series and expanding 

H (r) in a Bessel series leads to 
m 

b., = 2 s2 Cr)J. (S. .r)rdr/{(S.2 _ l2)j2(s )}. 

ij aj 0 1 a ij ij X xj 
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This Integral must be computed numerically. For the sake 
of brevity only the first five coefficients in each series 
are calculated and presented in Table 6. 

Tab le 6 

BESSEL SERIES COEFFICIENTS FOR ANALYICAL SOLUTION 


n 

^On 

^2n 

0 

1. 11738-0. 37246y 


1 

0.50462+0. 37912y 

0. 27181-1. 09482y 

2 

-0. 09246-0. 00802y 

-0. 11390-0. OIIIOy 

3 

0. 05019+0. 00184y 

0. 05479+0. 00208y 

4 

-0. 03290-0. 00067y 

-0. 03451-0. 00071Y 

5 

0. 02375+0. 00030y 

0. 02451+0. 00033y 


In Chapter 5, equation (5.8) can be written as 


^ ~ ^I'^’l ^ ^ 2^2 ^1*^3 (A-1) 

where 4 >q = Jq, (j)^ = J^cosQ, <t >2 = J 2 COS 20 

(}>2 = J^sin0, = J2Sin20 . 

By substituting (A-1) into the governing equation (5.6), 
yields 


R = D((f.) 


CA-2) 




where D is the nonlinear differential operator of equation 
(5-6). If (A-1) represented the exact solution to (5.6), R 
would vanish. Since this is not the case R has a finite 
value. The Galerkin procedure consists of making R orthog- 
onal to each of the <ti^’s. This leads to the equations 

/ V^’^R<|),rdrd =0 i = 0,1,. ..,4 (A-3) 

0 0 1 

and a set of five second order differential equations are 
generated. The angular part of integration can he performed 
in closed form while the radial part of integration must he 
computed numerically and leads to the integral 



<^1 = fk//ie)J^CS^jrIrdr 

(A-4) 

where 


(A-5) 


The numerical values of are listed in Table 7 for 
Y = 1.2. Also listed in this table are the functional forms 
of P^’s. In writing these the notation 

Rl = S2 ^(y- 1)JJ + - i2j^/r2, i=0,l,2 (A-6) 


is used. Similarly the numerical values of 



and functional-, forms of G 


ij 


s are listed in Table 8. 
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Table 7 


COEFFICIENTS APPEARING IN EQUATIONS (5-9) Y = 1-2 


i 

k 

^1 


1 

0 

4.137 

RqJq + 2ShJi^ 

2 

0 

1.042 

R^J^/2 + Sf^J]2+ jz/r^ 

3 

0 

-0.208 

R 2 J 2/2 + + 4j|/r2 

4 

1 

-1.939 

2 S 01 S 21 JJJI 

5 

1 

-2.312 

RiJq + 2SQ2S^2'^o'^i 

6 

1 

1.719 


7 

1 

1.483 

R 2 J 1/2 + + 2J^J2/r2 

8 

2 

-2.785 

R 0 J 2 + 2SqjS22Jq12 

9 

2 

-3.039 

Vo + ^^01^21^0^2 

10 

2 

1.132 

RjJj/2 + S2^J]2 „ j2/r2 

11 

0 

2.586 


12 

0 

0.480 

JsCShJi^ + ■Jf/r’2) 

1,3 

0 

-0.196 

221^2^/2 + 2j2/r2 

14 

1 

-4.849 

2S S J’J' 
01 11 01 

15 

1 

1.503 

S, S J'J’ + 2J J /r2 
11 21 1 2 12 

16 

2 . 

-0.576 

^^(S2iJj2 _ j2/p2) 

17 

2 

-3.481 

2S S J’J' 
01 21 0 2 
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Tab le 8 


COEFFICIENTS APPEARING IN EQUATION (5.12) 


■ 

1.000 

<1 0 CSoir) 


1.293 



0.240 

+ Jf/r2) 

- 

0.098 

"^(S2^J^2/2 + 2J2/r2) 

- 

0.176 

-hJl 


0.061 

1 

-VI 


0.049 

-V| 


1.000 

JiCSiir) 

- 

1.212 



0.927 

J5(S^^S2iJ|J> + 2J^J2/r2) 


0.165 

-JgJl 

- 

0.199 

2 J 2 


1.000 

J 2^321^) 

- 

1.741 

SoiS2iJ;J2 

— 

0.223 

Ji;CS2,J«2 _ j2/p2) 



APPENDIX B 


PROGRAM FOR STABILITY BOUNDARY CURVES 
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SS£T LIST FREE 
C 

C KIN-WING WONG JUN 30 1978 

C 

C CCMBUSTIQN INSTABILITY 

C 

C METHQO GALERKIN METHOD 

C« THIS PROGRAM GENERATE THE STABILITY CURVE. 
C 

C INPUT INFORMATION 

C 


c 

THE PROGRAM 

USED THE FORMAT FREE INPUT 

OUTPUT. 

A 

c 

COMMA SERVES 

AS A DELIMITER AND THE INPUT VARIABLES 

c 

AFE ENTER IN 

THE FOLLOWING ORDER 



c 

TI 

INITIAL TIME 



C 

TF 

FINAL TIME 



c 

DT 

STEP SIZE FOR TIME 



C 

WB 

BURNING RATE 



c 

p 

TOELAY 

DELAY TIME 



c 

El 

INITIAL EPS 



C 

EF 

FINAL EPS 



c 

DE 

STEP SIZE FOR EPS 



c 

DA 

ASSUME POSITION FOR N 



c 

FO»FD 

ARE THE INITIAL CONDITION 

MODIFIER 

AND 

c 


HAVING THE VALUE EITHER 1 

OR 0 



c 

C ANOTHER SETS OF DATA CAN BE ENTER BY ENTERING 

C EI» £F.>D£i»OA»FO AND FD . 

C 

DIMENSION DYaO)»FO(5),OF(5),G(10#500 ).Y<10) 

CCHMON EP* RK#hB>.EPANV 
COMMON /SLISr/S0^Sl»S2 

COMMON /BL«l/Cl#C2»C3»C4i.C5i.C6^C7#C8#C9i.C10»Cll»C12 
l.C13»C14. C15^C1S»C17 

LOGICAL PASS^GROh 
DATA KAXITR»EPS/20».l/ 

DATA PI/3.14159/ 

READ(5»/) TI»TF»DT»WB»TO£LAY 
FFIN T*//. TI^ TF#D T* wB.TDELAY 
TD£LAY=PI*TDELAY 
K=TOELAY/OT 

M = K41 

IF (K1.GT.5C0) GO TO 105 
If (K.LT.l.AND.NCF.EC.l ) K=1 
N STEP= CTF- TI )/D T 
N = 10 

HhALF=0T*.5 

NFR03=0 

9 REAO( 5#/»£ND = 99 9 ) E I . EF,. D£ »0A ^F0» CF 

NFR08=NPR08*1 
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FFINT*//»NPF0a 
PFINT *//#F0»DF 
MSTEP*(£F-ei )/0£ 

EF=EI 

IF (EP *LT.O) £P=.01 
FASS=,FALSE. 

CC IOC ISTEP =Os.MSTEP ' 

DC 80 ITR=1»MAXITR 
EFANV = £P*ANV 

IF (EPANV-EC-0) £PANV=ANV 
CC 10 I=1»K 
CC 10 J = l»N 
1C G(JrI)=0. 

CC 20 I=1»N 
2C V(I)=G. 

C 

C INITIAL CONDITION 

C 

y (2)=F0C1) lYCA ) = F0(2) J Y( 6)=F0 ( 3 ) I Y ( S ) ^F0< A ) I Y C 10) = F0 

1(5 ) 

Y(1>=DF(1)*S0;Y(3>=DF<2)*S1;YC5)=DF(3)*S2 
Y(7) = DF( 4)*S1;Y< S) = DF(5)*S2 
CALL DLYFUN ( N » Y » G EPS o Kl ) 

T = TI 
C 

CC AO I=i»NSTEP 

CALL BKDELY (N^T<.Yi»DY»OTpG»KpEP»8illOj>HHALf »«K1) 

CC 30 J=2pN»2 

1FCABS(Y(J )).GT.2) GO TO 60 
3C CCNTINLE 
AG CCNTINUE 

GFOW=. FALSE. 

If (PASS) GO TO 50 
AN YDWN = ANY 
ANV=ANV+OA 
GC TO 80 
50 CCNTINUE 

AAVOMN^ANU 
GC TO 70 
60 CCNTINUE 

PASS=. TRUE. 

GFOW=-TRU£. 

ANVUP=ANV 
7C CCNTINLE 

IF ( ABS(ANVLP-ANYDHN).LT-EPS) GO TO 90 
ANV= ( ANVUP + ANVD hN)*. 5 
ac CCNTINUE 
SO CCNTINUE 
NNU = AN V 

PRINT *//^EP»-NNU»ITR 
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IF ( ITR.EQ.HAXITR ) PRINT //»«*<** i^ARNING *** THE 
1 ANSWERMAY NOT CONVERGE* 

EF = EP-»DE 

ANV=ANVUP/2#ANVDWN=0 
IGC CONTINUE 
GC TO 9 

1C5 PFINT //**DELAY TIME TOO LARGE *»TDELAY 
GC TO 999 

lie FFINT *//»*ILLEGAL DELAY FUNCTION* 

9S5 STOP 
END 

SUBROUTINE DIFFUN (T»Y»DY) 

C 

C THS SECTION PROVIDE A SET OF FIVE SECOND ORDER 

C EQUATIONS 

C WHERE 

F0=Y<2) G1=Y<8) 

F1=Y<A) G2=Y(1C) 

F2=Y(6) 

DIMENSION Y{10)»CY(10) 

COMMON £P, RK»WB#EPANV 
COMMON /SLIST/S0/S1»S2 

COMMON /BLia/ C1»C2»C3^CA,C5»C6/C7 »CS»C9# ClCi.Cll»Cl2 
1»C13»C14» C15»C16#C17 

C 

CU2 )= Y(1 ) 

DY(1) = -S0*SC»Y(2)-EP*(C1«>Y(2)*YC1>*C2«(Y(4)»Y(3)^Y(8) 
1*Y(7)) ♦C3*(Y(10 )*Y(9) + Y<6)*Y(53>) 

CH4 ) = Y(3) 

DY(3)=-S1*S1*Y(A)-EP*'(CA*YC2)*Y(3)^C5»Y(A)*Y(1) ♦ 

1 C6*(Y(8)*Y(9)*Y<A)*Y(5 ) )+C7» < Y (6 )* Y ( 3 )♦ Y (1 0 )<» Y< 7 ) ) ) 
CY(6)=Y(5) 

CU5 ) = -S2* S2*Y(6 )-EP*(C8*Y(2)*Y(5)*C9*YC6 >*Y<1 )♦ CIO* 
1(Y(8)*Y(7)-Y(A)*Y(3))) 

CU8)^Y(7) 

CY(7) = -S1*S1*Y(8)-EP*(CA*Y(2)*Y(7) + C5*Y(8)*YU)4 C6* 

K Y( A) *Y(9)-Y(e)*Y(5))*C7*( YaC)*Y{ 3)-Y(6)*YC7 )) ) 

OYO 10)=Y(9) 

CU9)=-S2*S2*Y(10)-EP*(C8*Y{2)*Y(S)+C9*Y<10)*Y(1) CIO 
1*(YC 3)*Y( 3 )*Y(A )*Y(7n ) 

IF (RK.EQ. 2) RETURN 
C 

C THIS PROVIDE THE COMBUSTION TERMS 

C 

CY( 1)=DY(1)-WB*(Y(1 )+EPANV*(Cll*Y(2)*Y(2)4Cl2*{Y( A)*Y 
HA)* Y(8)*Y(8))* C13*0 YUO >* YUO )*Y<6 )*Y(6))>) 

CU3 ) = DY(3 )-w8*(Y( 3)*EPANV*(C1A*Y0 2>*Y0 A) + C15*-{ Y( A )*Y 
1(E)* Y(3 ) *Y( 1 0 ) ) ) ) 

DY(5) = DY(5)-Wa*(Y(5)*EPANV*(C16*{Y(A)*Y(A)-Y<8)*Y(8)) 

1 *C17*Y(2)*Y(6))) 
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CH 7) = DY(7 )-W8*(V(7)-»EPA NV*(C14*Y<2 )* Y(8 )^C15*< Y( A )*Y 
1(10) Y(8)*Y(6)))) 

QY( 9) = 0Y<9 )-Vi0*r(Y<9)+EPANV*<Cl?*Y<2)*Y(lO)+2*C16*Y{4> 

1 * Y ( 8 ) ) ) 

RETURN 

END 

StBRQlTINE RKDELY(N»Tf Y»DY/H»G>K»EPS» *^HHALF»Kl) 

C 

THIS SECTION PERFORM A FOURTH ORDER RUNGE*KUTTA METHOD 
NITH TIME DELAY FUNCTION 

DIMENSION Y(1Q) »DYn0)#Y2(10) »Y3t 10 )»G( 10^500) 

CCMNQN EP» RK^feB^EPANV 

CCMMQN /BLiU/ C 1 i»C2#C 3^ C 4 »C 5» C 6 j. C7 » C 8» C 9, Cl 0/ C 1 1 » C 1 2 
1»C13»C14» C15»C16^C17 

CALL DIFFUN(T#Y#DY> 

CG 10 I=1»N 

Y2 ( I )=Y (I ) ♦ HHALF*(DY( I > 4 G(I#D) 
G(I»l)=(G(I»l)4G(I»2))/2. 

1C CCNTINUE 
T=T4HHALF 

CALL CIFFUN(T»Y2#DY) 

CC 20 I = li.N 

Y3(I)=Y (I ) 4 HHALF*CDY<I ) ♦ G(I»D) 

2C Y2(I)= Y2U) * 2,^ Y3( I ) 

CALL OIFFUNC TpY3#DY) 

CC 40 I=1^N 

Y3(I)=YU) ♦ H*(DY(I) 4 G(I»D) 

IF(K.LT.l) GO TO 40 
CC 30 J=1^K 
3C G(I»J) = G(I » J4l) 

AC Y2(I) = Y2U) * Y3(I) 

T=T4HHALF 

CALL DIFFUN (T »Y3j>DY) 

CC 50 I=1*N 

YCI) = (Y2(I) - Yd) 4 HHALF*(DY<I) 4 6<l4l)))/3- 
5C CONTINUE 

ENTRY DLYFLN ( N » Y » 6» K, EP S»K1 ) 

C 

C THS GIVES THE FORM OF DELAY FUNCTION G 

C 

IF (K.LT.l) RETURN 
IF (RN.EQ.2) RETURN 

G(1»K1)= W8* EPANV*<C11*Y(2)*Y(2)4C124(Y(4 )*Y 

ld)4 Y(8)*Y(8))4 Cl3*(YaC)*Y(10)4Y(6)*Y(6))) 

G(3»K1)= we* EPANV*<C14* Y(2 )* y (4 )4C15*( Y( 4 )*Y 

1(6)4 Y<8)*Y(10))) 

G(5*tU) = WB*EPANV*<Cl6*(Y(4)*Y(4)-y(a)*Y<8))4C17*Y(2)*Y 
1(6 ) ) 

G(7»K1)= W8* EPANV*(C14*Y(2)»Y(8)4C15*C Y(4)*Y 

1(10) Y(8)*Y(6))) 
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G(9»K1 )=W8*EPANV*(C17*y( 2)*Y( 10)-»2-kC16»Y(4>*Y(e)) 

SETURN 

END 

BLOCK DATA 

COMMON /8LK1/C C2 » C3 »C A »C 5 *0 6»C7 » C8» C9» C 10 »C 11 »C 1 2 
1#C1 3#C14» C15»C16#C17 

CCMMON /SLIST/S0»S1^S2 

DATA S0»S1 »S2/3 -83 171 >1.8 4 118 >3.05424/ 

CATA Cl>C2>C3>C4/4,1373>l.0423>-.2084>-1.9394/ 

DATA C5>C6>C7>C6>C9 /*=2. 312 3> 1-7187 >1 .4828>-2- 785> 
1-3.0386/ 

CAT A C10>C 11>C12»C13/1.1318>2 .586>.480>-. 196/ 

CAT A C 14>C 1 5>C 1 6>C 17/-2. 424 3> 1-353 4 44 5>-. 447 >-3.481/ 
END 



APPENDIX C 

PROGRAM FOR MODAL AMPLITUDE AND PRESSURE 

■) 
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SSET LISTUTQ8IN0 

S3IN0 = FROM *F0RTLI8/= FREE 

FILE 1(KINC=DISK»MAXRECSIZE = 15»8L0CKSIZE=A20» AREAS=<» 
AREASI2E = 1A00»FILETYPE = 7 ) 

C 

C KIN-WING WONG JUN.30 197Q 

C 

C CCMBUSTION INSTABILITY 

C 

c method; this program lsed the galerkin method to 
C generate FIVE SECOND ORDER DIFFERENTIAL EQUATIONS. 

C A FOURTH ORDER RUNGE-KUTTA METHOD IS USED TO SOLVE 

C THE RESULTING EQUATIONS. 

C 

C THE PROGRAM GIVES THE MODAL AMPLITUDE AND THE PERTUBED 

C PRESSURE IN GRAPHICAL FORM. 

C 

C THERE ARE TWO TYPES OF INPUT DATAS. THE FIRST SET OF 

C INPUT data used FORMAT FREE INPUT. A COMMA SERVES 

C AS A DELIMITER AND THE DATA ENTER AS FOLLOrt; 

C TI INITIAL TIME 

C TF FINAL TIME 

C OT STEP SIZE FOR TIME 

NP PRINT FREQUENCY 

NDF 1 FOP THE DELAY TIME APPROCHING 0 
0 OTHERWISE 

IGAS I FOR GAS-DYNAMIC TERMS 
0 OTHERWISE 

F0»' INITIAL CONDITION MODIFIER 
OF El THER 1 OR 0 
R RADIUS OF THE CYLINDER 

THE SECOND SET OF DATA USED THE NAMELIST (LIST) 

INPUT OUTPUT OPTION. THE DATA CAN 8£ ENTERED IN 
ANY order. however THE FOLLOWING DATA MUST 8E 
PROVIDED INITIALLY. 

TDELAY DELAY TIME IN MULTIPLE OF PI 
EP ORDER PARMETEH 
W8 STEADY BURNING RATE 
ANV INTERACTION PARMETER 
PLOT T FOR PRESSURE PLOT 
F OTHERWISE 

IF PLOT IS TRUE» ONE MORE DATA (NUMPT) IS NEEDED TO 
PROVIDE THE NUMBER OF POINT WITHIN 0 AND 2*PI. 

C 

DIMENSION XX(lO0)»YY(lOO),OY(lO)»P(5)^FO(5)»DF(5)»G 
1(10 j-5C0)»Y(10) 

CCMMON EP»IGAS/ RK, W3.EPANV 
CCMMQN /SLIST/S0,S1»S2 

CCMMON /BLKl/ C UC 2 ^C 3 » C A ^ C 5i- C6 ^C7 » C3 ^C9 » Cl 0^ Cl 1 » C 1 2 
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1»C13,CIA» C15#C16*C17 

COMMON /BLK2/C02;>C03»C04»C05^C06#C07,C12P#C1 3P»C14P 
1>C15P» C22^C23»C2-!.»C25,GAMA 

COMMON /8LK3/ 3 J 0» B J1 » B J2 » COS Q»C QS2 Q» SINQ# S IN 20 
DATA PI/3-1A159/ 

LCGICAL PLOT 

NAMELIST /LIST/ TDELAYi>EP»W8> ANV^PLOT 
READ (5,/) TI »TF» DT#NP»NOF»IGA S#FO;.OF»R 
PRINT *//#TI#TF»DT»NP»NDF 

IF (IGAS.EQ.O) PRINT //»«GAS DYNAMIC TERM IS QFF» 

PRINT *//»FO»OF»R 

EFS=l.E-5 

CALL BESJ( SO*R# 0^8J0^EPS*IER) 

CALL 3ESJ( S1*R» 1»BJ1»EPS»1ER) 

CALL BESJC S2*R^2»8J2^EPS»IER) 

999 R£A0(5#LIST»END=99) 

REWIND 1 

TCLY=PI*TDELAY 

EFANV=EP*ANV 

IF (EPANV.EQ.O) EPANV=ANV 
WFITE(6»U ST) 

N = 10 

riHALF=DT*.5 

K=TOLY/DT 

K1=K-K 

IF (K1 ,GT. 500) GO TO 500 
IF (K.LT,l.AND.NDF.EQ.l) K=1 
DC 30 I=1^K 
DC 30 J=1»N 
30 G(J»I)=0. 

DO 10 I=1?N 
10 Y(I)=0. 

C INITIAL CONDITION 

C 

T=TI 

Y(2)=F0(1))Y(A)=F0(2)JY(6)=F0(3); Y(3)=F0(4);ya0) = F0 
1(5) 

Y(1)=0F(1)*S0;Y( 3)=DF(2)*SI;Y(5)=CF(3)*S2 
Y(7)=DF(A) *S1;Y{9)=DF(5)*S2 
NSTEP=(TF-TI )/DT 
IF (.NOT. PLOT) GO TO AO 
CALL PRE3UE(Y^P) 

WRITECI) T»P 
AC CONTINUE 

PRINT //^"FUNCTIQN" 

WRITE(6.100) T^( Y(I ).I = 2,-N#2) 

A5 CONTINUE • 

CALL OLYFUN ( N » Y »G » K » EPS » K 1 ) 

DC 20 I=1^NSTEP 

CALL RKO£LY(N»T» Y»OY.OT»G»K»EP> &A00>HHALF»K1 ) 

DC 55 J=2»N»2 
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IF<A8S(Y(J)).LT.2) GO TO 55 
R^ = 2 

55 CONTINUE 

IF (MQO(I^NP).NE.O) GO TO 20 
IF (.NOT.PLOT) GO TO 60 
CALL PRESU£(Y»P) 

WFITECi) T»P 

60 xFITE(6#100) T#(Y(J)#J=2#N»2) 

65 CONTINUE 
2C CCNTINliE 

IF (PK.GT.O) print // » » UNST A3LE* 

IF (, NOT, PLOT) GO TO 70 
ENDFILE 1 
REWIND 1 

PRINT //»» PRESSURE * 

REA0(5,/) NUMPT 
OS=PI *2./NUMPT 
NLM=NUMPTt 1 

TEM°=c;xx{ i)=o; 

DC 75 I=2»NUM 
T£MP = T£MP-» DO 
75 XX(I)=TEMP 

98 R£A0(l#EN0=70) T^P 
PRINT //»"TIM£= "#T 
SUMY=C. 

DC 80 1=1# NUM 

ae=xx(i ) 

YVT I )= “GAMA*{ P{ 1)<^P(2) »C0SCQQ)4p{3 )*COS (2»QQ)4P<4 ) 
1*SIN(QQ)4 P(5)*S1N(2*QQ) ) 

SUMY = SUMY» YY(I) 

80 CONTINUE 

IF (SUMY.NE.O.) CALL PLOT 20 ( X X # Y Y#NUM # 1 00 # 1 00) 

GO TO 98 
70 CCNTINUE 
RK=0 

GC TO 999 

400 PRINT *//,»ILL£GAL DELAY FUNCTION* 

GC TO 99 

50C PRINT *//»*0£LAY TIME TOO LARGE* 

99 STOP 

100 FORMAT (X#8E12.5) 
lie F0PMAT(13X#7E12.5) 

END 

SU8R0UTINE OIFFUN (T#Yi-OY) 

DIMENSION Y(10)#OY(10) 

COMMON £P#IGAS# RK#W3#EPANV 
COMMON /SL IST/S0#S1#S2 

COMMON /8LK1/ C 1 #C 2#C 3 # C 4# C 5# C 6 # C7 # C 8 #C 9 # Cl 0 »C 1 1 > C 12 
1»C13#C14# C15#C16#C17 

C 

CY(2)=Y(1) 
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D Y<i )=-S0*S0*Y(2 )-£P<f IGAS*(C1*Y<2)*Y( l)*C2*(Y(4)ftY(3) 
1*Y(3)*Y<7)) »C3*(Y(10)*Y(9)4YC6)*Y(5))) 

DY(4) = Y(3) 

DY(3 ) = -Sl*Sl*Y(A )-EP*IGAS*(C4*Y(2)»Y( 3)*C5*Y(4)*Y<1 ) * 

1 C6*(Y(8)*Yt9)+Y(4)»Y(5))4C7*(Y<6)*Y(3)4YC10)*Y 

2 £ 7 ) ) ) 

0Y(6)=Y(5) 

DY(5 ) = -S2*S2*Y( 6 ) “E P* 1 GA S * ( C8 * Y ( 2) * Y( 5 ) 4 C 9*Y ( 6 ) *Y( !)♦ 

I C10*( Y(8) *Y(7 )-Y(4)*Y(3 )) ) 

DY(6)^Y(7) 

DY(7) = -S1*S1*Y<8)-EP*IGAS*(C4*Y(2)*Y(7)4C5*Y(3)*Y(1)* 

1 C6*{Y(4)*Y£9)-Y£8)*Y(5)) + C7*(Y£10)*Yt3)-Y£6)*Y<7))) 
DY( 10 ) = Y(9 ) 

OYC 9 )=-S2*b24Y( 10)-EP*IGAS*(C8*Y(2)*Y<9)4C9*Y(10)*Y£1) 
1 C10»(Y£8 )*Y(3)4Y(4)»Y(7))) 

IF (RK.EQ. 2) RETURN 

DY<1)=DY(1)-W8*(Y<1)*EPANY*(C11*Y(2)*Y<2)*C12*{Y(4)*Y 
1CA)4 Y(8)*Y(8)>4 C13*<Y(10>4Ya0)*Y(6 )*Y(6)))) 

CY< 3 ) = DY( 3 )-W8«r( Y£ 3)4EPANY *(C 14*Y£2) »Yt 4) + C15*( Y( 4 )*Y 
1(6)+ Y(8)*Y(10)))) 

CY(5 )=0Y(5 )-W8* ( Y(5)+£PANV*(C 16*£ Y( 4) * Y ( 4 ) -y ( 8 ) *Y ( 8) ) 

1 +C17+Y(2 ) + Y<6))) 

DY( 7)=0Y(7 )-wB*( Y(7 )+EPANY*(C14* Y(2)+y{8)+C15*(Y(4)+Y 
1(10) Y(8)+Y(6)))) 

CY(9 ) = 0Y(9 )-W8* ( Y(9 )+EPANV*(C17*Y(2)*YU0)+2 + C16*Y( 4) 
1*Y(8 ) ) ) 

RETURN 

END 

SU3R0UTINE R KD EL Y<N # T» Y /• 0 Y ? G *> Ko£ PS» * p HHALF *-K 1 ) 
DIMENSION Y(10 )^DY(10)»Y2(10) »Y3( 10 >»G( 10 » 500) 

COMMON EP»IGAS» RK»W8»EPANV 

CCMMON /0LK1/ C1»C2»C3»C4^C5^C6>C7»C8»C9»C10,C11»C12 
1»C13^C14» C15»C16».C17 

CALL DIFFUN(T»Y^DY) 

OC 1 1=1, N 

Y2(I) = Y(I) + HHALF*(DYU) + G(I,1)) 

G(I,1)=(G( I,l)+G(I,2))/2 

1 CONTINUE 

T = T+ HHALF 

CALL OIFFUN(T, Y2,DY) 

DO 2 1=1, N 

Y3(I)=Y(I) + HHALF*(DY(I) ♦ G(I,D) 

2 Y2(I)= Y2(I) ♦ 2*Y3(I) 

CALL DIFFUN( T,Y3»DY) 

DC 3 1=1, N 

Y3(I) = Y(I) + H + (DY(I) ♦ G(I,D) 

IF(K.LT.l) GO TO 3 
DC 10 J = 1,K 
10 G(I, J ) = G(I , J+1 ) 

3 Y2(I)=Y2U) + Y3(I) 

T=T+HHALF 
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CALL DIFFUN CT ,Y3»DY) 

DC 4 1=1»N 

Yd) = (Y2(I) - Yd) ♦ HHALF*(DY(I) ♦ Gd>l)))/3. 

A CONTINUE 

ENTRY OLYFUN ( N » Y G> K» EPS > K1 ) 

IF (K.LT.l) RETURN 
IF (RK.EQ.2) RETURN 

G<l»Ki)= W3* EPANV*<C 11*Y(2)*Y(2)*C12*( Y( A)*Y 

UA)+ Y(8)*Y(8))+ C13*(Y(10)*Y(10)4Y(6)*Y(6))) 

G(3»K1)= W8* EPANtf«>(ClA*YC2)»y<4)4C15*CY(4)*Y 

1(6)4 Y(8)4YC10))) 

G(5»K1 ) = W8*EPAN V*(C16* (Y( A)*Y(4)-Y( 8) *Y ( 8 ) )♦ C 1 7* Y ( 2 ) * Y 
1(6) ) 

G(7fKl)= W8* EPANV*(C14*Y(2)*Y(8)4CI5*( Y(4)*Y 

1(10) Y(3)*Y(&>)) 

G(9»K1 )=Vi8*EPANV*(Cl 7*Y(2)*Y( 10) + 2*C16*Y( 4)*y<8 )) 

RETURN 

END 

SUBROUTINE PRESU£(Y»P) 

C 

C THIS SECTION PERFORM THE PRESSURE CALCULATION 

C 

DIMENSION Y(10)»P(5) 

COMMON EP^IGAS^ RK^WBj-EPANV 

CCMMON /8LK2/C02j.C03#C04#C05,C063>C07»Cl2P»Ci3P»C14P 
1»C15P» C22»C23^C24 #C25»GAM A 

COMMON /8L.-(3/BJ0<.aJl»8J2»C0Sa »C0S2Q?SINQ»-SIN2a 
OFO = Y( D; DF1 = Y( 3 )^ 0F2= Y ( 5 ) I DG 1= Y ( 7 ) 7 3 6 2= Y ( 9 ) 

F0^Y(2)> F1=Y(4)I F2 = Y(6)» Gl = Y{a>; G2=YdO) 

P(1 ) = e J0*( DF04EP*(C02*FO*FO-»C03*(Fi»Fl4Gl*Gi )4C04* (F2 
1*F24G2*G2) 4C05*OFO*DF04C06*{DF1»QF14DG1*OG1 )4C07* 

2(CF2*DF24 DG24DG2))) 

P(2)=3J14(0F14EP*(C12P*F0*F14C13P*(F1»F2461*G2)4C14P 
1*0F0*DF14 C15P*(DFl40F24DGl*DG2)> ) 

P(3)=BJ2*(DF24EP*(C22*F0 *F2 4C2 3 » (F 1 *F 1 - G1 *G1 )4C24*DF0 
l*0F2+ C25*( OF1*OF1-OG1*OG1 ) ) ) 

P(4 ) = B Jl*( DG14£P*(C 12P*F0*Gl4C13P*(Fl *G2-G1*F2 >4C14P 
1*CF0 ‘■DG14 C15P*(0F1*0G2-DF24D61 ) ) ) 

P (5 ) = BJ2*{ OG2+EP*(C22*FO*G242.*C23*F1 *G1* C24»0F0*DG24 
1 2-*C25*DF1*DG1 )) 

RETURN 

END 

3L0CK DATA 

COMMON /BLKl/ C 1 .-C 2»C 3, C 4> C 5^ C 6# C7 r C8 ^C 9 » Cl 0 1 1 ^ C 12 

1»C13»C14» C15/C16»C17 

COMMON /SL I ST/SO irSl»S2 

CCMMON /BLK2/C0 2»C0 3rC04rC0 5rC0&j-C0 7»C12P*C13P»C14P 
1»C15P» C22^C23i-C24 fC25»G AM A 

DATA SO^Sl » S2/ 3. 831 71 1.841 13 » 3. 05424/ 

DATA C1»C2#C3»C4/4,1373» 1.0 42 3j»-.20 84^-1.9394/ 

DATA 05^06 rC7»C8»C9 / - 2. 31 2 3» I . 7 187 »1 . 4 828 » -2 . 785 » 
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1-3.0368/ 

DATA C iO.C 1 1»C12>C13/1.1318.2 .586». A80»-. 196/ 
DATA CIA»C15»C16#C1 7/- 2 . 42 A 3 r 1.85 3A 4 A 5 ^ -- UAl f 
DATA C02»C03#C0A»C05»C06»C07/ 1.^930»0-^AOO^-- 
l-.1762 >.0607,.0A9A/C12P»C l3Pj»Cl AP .C15P/-1 

2» .9267 ,.1651 1907/C22»C23»C 2A»C25/- 1.74 06»- 
3, .2371#-. 1754/ 

DATA GAMA/1.2/ 

END 


-3-481/ 

0982# 

.2121 

.2235 
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